Attraction Between Like-Charged Walls: 
Short-Ranged Simulations Using Local Molecular Field Theory 
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Effective attraction between iike-charged waifs mediated by counterions is studied using iocaf 
mofecular fiefd (LMF) theory. Monte Carfo simuiations of the "mimic system" given by LMF theory, 
with short-ranged "Couiomb core" interactions in an effective singie particle potential incorporating 
a mean-field average of the long-ranged Coulomb interactions, provide a direct test of the theory, and 
are in exceilent agreement with more complex simulations of the full Coulomb system by Moreira 
and Netz [Eur. Phys. J. E 8, 33 (2002)]. A simple, generally-applicable criterion to determine the 
consistency parameter a m i n needed for accurate use of the LMF theory is presented. 



Effective attractions between like-charged objects are 
uite common and have been extensively studied 0,0,0, 
For example, highly charged DNA is densely packed in 
cell nuclei via positively-charged intermediaries, and, in 
vitro, DNA may be condensed into toroids by adding suf- 
ficient concentrations of divalent or trivalent counterions. 
In this paper we use local molecular field (LMF) theory 
to study one of the simplest models that exhibits like- 
charged attraction — two uniformly-charged walls with 
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(a) Full Model System 



neutralizing point counterions, as shown in Fig. 1(a) with 
length scales that will be discussed later. 

LMF theory defines a general mapping that relates the 
structure and thermodynamics of a nonuniform system 
with long-ranged intermolecular interactions in an exter- 
nal field 4> to those of a simpler "mimic system" with 
short-ranged interactions in the presence of an effective 
field 4>r, as qualitatively shown in Fig. 1(b) 4>r accounts 
for the averaged effects of the long-ranged interactions 
and self-consistently depends on the nonuniform density 
the field induces 0, This approach is particularly 
useful for systems with Coulomb interactions, because 
one can choose specific slowly-varying, long-ranged com- 
ponents of the Coulomb interactions that are especially 
well-suited for the mean-field average. The remaining 
short-ranged "Coulomb core" components combine with 
other existing short-ranged interactions to define the in- 
termolecular interactions in the Coulomb mimic system. 
Thus the theory is not restricted to point counterions and 
very accurate results have already been found for uniform 
fluids with charged hard cores [3, la] • 

However, additional approximations were made in 
these earlier applications of LMF theory. In particular 
the Boltzmann approximation for the density response 
to the effective field was used for the charged wall sys- 
tem [ij]. Here we use Monte Carlo (MC) simulations to 
accurately determine the density response. We believe 
that such simulations of the mimic system will often be 
needed to obtain quantitative results from LMF theory 
in more realistic models of biophysical interest. 
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(b) LMF Mimic System 



FIG. 1: Ion distributions and length scales for a moderate 
coupfing case where the attractions are just beginning to de- 
velop (£ = 20 and d = 10). The full model system shown 
in |(a)| consists of point counterions represented by points that 
neutraiize two charged hard walls, with dashed circles for Lb 
(only 3 Lb circles are shown for clarity) and dotted lines for 
Lg- The electrostatic potenti al (f> due to the walls is 0. The 
LMF mimic system shown in |(b)| has Coulomb core interac- 
tions with a range of o m in proportional to the spacing Lw 
(indicated by solid circles), and a modified wall potential <f>a 
that accounts for the remaining long-ranged interactions. 



Then the only remaining errors are those inherent in 
the LMF mapping itself. The results provide a criti- 
cal test of LMF theory in a nonuniform Coulomb sys- 
tem where the basic physics of the counterion-mediated 
attraction is highly nontrivial but well-understood, and 
where extensive benchmark simulations are available HJ. 

Before giving details of the model and the LMF map- 
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FIG. 2: Pressure vs. distance curves for two couplings, £ = 
20 and £ = 100. LMF simulations agree very well with full 
simulations from 0. For d = 20 and £ = 20, the effective 
= 18 and for d = 20 and $ = 100, cr min = 34. 
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ping, the quantitative accuracy achieved in practice is 
illustrated in Fig. [21 which compares the dimensionless 
osmotic pressure P for the full Coulomb system |9( to 
results of the LMF theory for two system couplings, £. 
Depending on the separation d, either coupling can re- 
sult in a net attractive force on the walls as indicated by 
negative values of P. Simulations of the full Coulomb 
system required careful and costly treatment of periodic 
boundary conditions using the Lekner-Sperb method; our 
simulations of the short-ranged mimic system used only 
a simple minimum image method. 

In the model two infinite hard walls with a negative 
charge density q w are located at z = and z = d. Pos- 
itive neutralizing point counterions with valence Z and 
charge Ze$ are contained in < z < d and there is a uni- 
form dielectric constant e everywhere. The three length 



scales shown in Fig. 1(a) can be motivated by an exami- 
nation of the energetics of one uniformly charged wall at 
z w . The potential energy between the wall and a counte- 
rion at z is —2Trq w Zeo \z — z w \ /e. The distance where 
this potential equals ksT defines the Gouy- Chapman 
length (including ion valence) Lq = efceT l /\2irq w Zeo\. 
The Bjerrum length L b including ion valence is similarly 
defined using the potential energy between a pair of ions: 
Lb = Z 2 el/(eksT). The third length scale, L w , is de- 
termined from the surface area of the wall neutralized by 
one counterion: L? w = Ze / \q w \ — 2itLbLq. 

We will use dimensionless variables where lengths are 
measured in units of Lq and energy in units of UbT. 
Specifying d and the coupling strength £ = Lb/Lq fully 
defines the thermodynamic state of this system. Effective 
attractions can arise for strong-coupling states with £ > 
12 0. Since the total force on a counterion from both 
walls exactly cancels, the bare external potential <j> — 
for < z < d. However due to the long-ranged Coulomb 



repulsion, counterions will organize next to the walls into 
either one or two layers, based on a complex balance 
between coupling strength £ and the width d available. 



Figure 1(a) qualitatively depicts a weakly attractive 
state with £ = 20 and d = 10. Here, most counterions are 
found in separate two-dimensional (2D) liquid layers near 
each wall, with a characteristic nearest neighbor spacing 
of order L w = 11.2 fixed by local neutrality. The effec- 
tive attractions arise mainly from cross-correlations be- 
tween ions in the two layers, and become even stronger at 
smaller separations when the counterions are forced into 
a single 2D layer with characteristic spacing L w j\/2 



Figure 1(b) gives the LMF mapping to the mimic sys- 
tem. Mimic ions have a short-ranged repulsive Coulomb 
core with size a min of order L w (indicated by solid 
circles), and their averaged long-ranged repulsions lead 
to an effective external field 4>r{z) with wells of depth 
l.hksT at the walls. At still larger d, the wells deepen 
and mimic particles form two distinct layers. At smaller 
d, the wells disappear and the Coulomb cores from sep- 
arate layers overlap significantly, forcing the system into 
a single 2D layer. 

The derivation of LMF theory and its application 
to general Coulombic systems are given in detail else- 
where |E 0, Si- The basic ideas can most easily be 
seen for a simple one-component system in an exter- 
nal field </>(r), where the intermolecular pair potential 
w(r) — uo(r) + u\{r) is properly separated into a short- 
ranged "core" part u Q (r) and a slowly- varying longer- 
ranged U\{r). The mimic system is composed of pair 
interactions Uq{t) and a renormalized or effective exter- 
nal field (pn(r), which is supposed to induce a nonuniform 
singlet density in the mimic system (indicated by the sub- 
script R) equal to that induced by <f> in the full system: 



Pfl(r; [<Pr]) = p(r; [4>] 



(1) 



This defines a mapping relating structure in the mimic 
and full systems. The explicit LMF equation for </>_r(i") 
incorporates a density-weighted average over the slowly- 
varying interactions Ux and is given up to a constant by 



<M r ) = fa 



drVfl(r';[fo])«i(|r'-r|). (2) 



This equation can be derived by integrating the first 
equation of the exact Yvon-Born-Green hierarchy relat- 
ing intermolecular forces and conditional singlet densities 
in the full and mimic systems after making two intercon- 
nected and physically reasonable approximations 0, ||| . 
First, when Eq. (JTJ holds, the conditional singlet density 
pn(r'\r; [</>_r]) in the mimic system should also approxi- 
mately equal that in the full system, provided that Uq 
gives a good representation of the short-ranged core in- 
teractions between particles. Second, the force from the 
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slowly-varying ui(|r' — r|) should be very small over the 
range a of characteristic nearest neighbor distances where 
p R {r'\r; [<j> R ]) differs significantly from p R (r'; {<f) R }). Then 
p R {r'\r; [<j> R ]) may be reasonably replaced by p R (r'; [</> R ]) 
in the integration of the force that yields Eq. J2J. 

For Coulomb systems we can control the accuracy 
of the second (mean field) approximation by convolut- 
ing the dimensionless Coulomb potential w(r) = £/r 
with a Gaussian whose width a is a parameter at our 
disposal 0, 01 ■ This yields a long-ranged component 
iti (r) = £erf(r/<7)/r that remains slowly varying at dis- 
tances less than a as illustrated in Fig. 3(a) and the as- 
sociated core component uo(r) = £ erfc(r/cr)/r. When a 
is too small, results of the LMF theory will be poor and 
will vary rapidly as a increases. But when a > a m im 
with drain of order the characteristic neighbor spacing 
a, we expect that u\ is sufficiently slowly- varying that 
the LMF averaging is consistent and there will be little 
change in results as a increases beyond <J m in 0- • 

Using this choice of u\, and noting that cf> = 0, we can 
integrate exactly over lateral coordinates in Eq. J5J) and 
obtain the two-wall LMF equation: 



dz'n R (z';[<t> R ])G(z',z). 



(3) 



Here n R (z) = 2n^p R (z) is a dimensionless rescaled 
density and G(z',z) = —\z — z'\eii (\z — z'\/ a) — 
(T7T -1 / 2 exp [— (z — z') 2 /a 2 ] + C(z') can be interpreted 
as the potential at z due to a Gaussian charge density 
CT7T -1 / 2 exp [— (z — z') 2 /cr 2 ] , with C(z') chosen so that 
G(z', 0) = 0. As explained earlier, (j> R (z) plays an im- 
portant role in this nonuniform mimic system. The bare 
field 4>{z) = must be replaced by 4> R {z) in simulations 
of the mimic system for the particles to separate correctly 



into two layers, as shown in Fig. 3(b) 

LMF theory gives exact results both as £ — > 0, where 
it reduces to the Poisson-Boltzmann (PB) theory, and 
in the strong coupling limit £ — > oo To assess its 
performance for intermediate £, canonical MC simula- 
tions of the nonuniform mimic system were carried out 
in a simulation cell of volume L x L x d. Particle num- 
ber N was chosen so that the wall length L dictated 
by neutrality is large enough to justify the minimum 
image convention. We used the simplest self-consistent 
simulation closure of the LMF equation by explicitly 
iterating solutions indexed by i of Eq. @ with con- 
verged MC values for n R (z) until the self-consistency cri- 
terion J dz 4> R — /d < 0.001 was met. Properties 
were averaged over 5 x 10 5 - 2 x 10 6 simulation sweeps. 
n R (z) and (j}Ri z ) were calculated on a grid with spac- 
ing Az = min {0.1, 0.Old}. The reduced osmotic pres- 
sure was calculated by an accurate method that uses the 
well-converged density at the midplane and the force be- 
tween particles and walls on the left and right of the 
midplane H M EH : P = n (z mldplane ) + 2^ (F LR ) /A. 



=£ 1 









lr 








\ \ 




— u (r) 
.... Ul (r) 




U.o 




/ \ 








- w(r) 




0.6 
0.4 








\ \ 












\ 








— n (z;[« 




\ 


a 




0.2 




.■■.n R (z; WR ]) 


V r 






- 4. _ 









10 20 30 

r 



(a) Pair Potential Split 



5 10 15 

z 

(b) Impact of <pji 



20 



FI G. 3 : Potentials and densities for £ = 10 wit h a = 14. 
In |(a)| £/r is split into uo and u\. As shown in |(b)| for d = 20, 
when simulating mimic particles that interact only via uo , the 
inclusion of (j>n, which here has well depths greater than 5ksT 
(scale on right) rather than the flat cj>, is crucial. no(z; [<f>]) is 
quite different than the accurate nn(z; [4>r]) as shown by the 
contact densities (scale on left) . 



Crucial for the success of the theory is the proper 
choice of <J m i n , which should scale with the characteristic 
neighbor distance a at strong coupling 0, Q • The dis- 
cussion above suggests a simple criterion that uses the 
consistency of the theory itself to determine a precise 
value for a m in- During a simulation with a given er, we 
measure the nearest-neigbor distance averaged over par- 
ticles, (L nn ). As a increases from small values by steps 
j, we expect initially that (L nn ) will increase as the core 
repulsions in uq(t) becomes larger. But for a > a m i n , the 
variation in (L nn ) should level off. As a numerical crite- 
rion that gives very reasonable results we choose the first 
a such that {{U nn ) ~ (Li- 1 )) / {a* - a^ 1 ) < 0.005. 

Figure 01 illustrates the application of this convergence 
criterion to a strongly- coupled system with £ = 100. In 
Fig. 4(a) the smoothed curves qualitatively represent the 
limits of amin determined for two separate layers at large 
d and for a single layer at small d. Between those curves 
lie data for specific d, ranging from two weakly corre- 
lated separate layers with weak attraction (d = 20) to 
a single layer configuration where attraction is maximal 
(d = 6). As the walls are pushed closer together and 
the counterions shift from two layers to one layer, the 
characteristic neighbor spacing a decreases by a factor of 
\/2. In the a m i„ convergence plots, we see the expected 
corresponding factor of y/2 as cr m in shifts from 36 to 26. 
shows how 



Figure 4(b) 



(L nn ) and a m i n relate to a 
special correlation function proposed by Rouzina and 
Bloomfield to better probe the nature of correlations be- 
tween particles in the 2D layers Q. The <72_d(?"||) shown 
here is the correlation function for pairs of particles on 
the same side of the midplane, with distances projected 
onto the a;y-plane. Note that a m m is larger than the 
distance a of the first peak in <?2d(?"||), ensuring that u\ 
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(a) a Convergence 



(b) Comparison to 
Correlations 



FIG. 4: Determination of the consistency paramet er a min . 
The criterion is applied for £ = 100 and varying d in |(a)| the 
horizontal line indicates the convergence threshold of 0.005. 
The smooth curves schematically represent the upper and 
lower limits for a convergence; once a single layer or two lay- 
ers unambiguously form, the neighbor distance a and hence 
drain should not vary further. (L nn ) (dotted lines) and ( Tjnm 
(dashed lines) are related to g2D pair correlation curves in|(b) 
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FIG. 5: Points where P(£, d) — from LMF simulations are 
compared to full system simulations in 0. In addition, sep- 
arations dmin where the maximum attractive force is found 
are fit by a power law to Lw oc (dmin = 1.18Lw 509 ). This 
agrees very well with an LMF scaling prediction 



varies slowly over this nearest neighbor distance. 

The "phase diagram" of points where P = 

is given in Fig. [SJ LMF results again agree very well 
with those of the full system. Also shown are points 
where LMF simulations yield the minimum pressure P m in 
(maximum attractive force). As argued in Ref. 0, P m in 
should occur at a separation d m in where a single layer of 
counterions first forms. By making a simple approxima- 
tion for (f)R and physically connecting a m i n to a oc L w , 
we find d m i n ck y/L w . The simulation results for d m in 
are best fit by d m i n — 1.18L?£ 09 , which agrees quite well 
with this scaling prediction [8j. 

Neither the use of an effective field in LMF theory nor 
simulations with short-ranged Coulomb cores is new. Im- 
plicit in the classical PB or Gouy-Chapman theory is an 
effective field resulting from a mean-field average of the 
Coulomb interactions. However, the PB approach in- 
cludes the rapidly varying Coulomb core components in 
the average (effectively choosing a = 0) and is accurate 
only for dilute weakly-coupled systems. Spherical trun- 
cations of Coulomb interactions as suggested by general- 
ized reaction field methods have had some notable suc- 
cesses in simulations of dense strongly-coupled uniform 
systems |l2j . But truncated interactions alone give poor 
results for geometrically nonuniform systems like water 
between wall s flii, Il4j or our two-wall system, as illus- 
trated in Fig. |3(b)| 

In contrast, LMF theory provides a general conceptual 
framework for nonuniform Coulomb systems. It deter- 
mines a physical choice for the short-ranged Coulomb 
cores and uses mean field theory in a consistent way to 
generate an effective potential that accounts for the re- 
maining long-ranged interactions. The classical PB ap- 



proach is greatly improved by averaging only over the 
slowly-varying u\ as dictated by LMF theory pfl and 
can even predict attraction for the two- wall system ||. 
We have shown here that simulations of the short-ranged 
cores in conjunction with the <pR given by LMF theory 
give quantitative agreement with simulations of the full 
Coulomb system. Further results for this system includ- 
ing analysis of the 2D correlation functions during the 
formation of a single or two layers, the scaling of a m in, 
and more realistic descriptions of counterions and co-ions 
will be reported elsewhere. This work was supported by 
NSF through grant CHE05-17818. JMR was supported 
by an NDSEG fellowship. 
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